Energy Barriers and Activated Dynamics in a Supercooled Lennard-Jones Liquid 
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We study the relation of the potential energy landscape (PEL) topography to relaxation dynamics 
of a small model glass former of Lennard-Jones type. The mechanism under investigation is the hop- 
ping betweem superstructures of PEL mimima, called metabasins (MB). From the mean durations 
(r) of visits to MBs, we derive effective depths of these objects by the relation i5 app = din (r) /d/3, 
where [3 — 1/feT. Since the apparent activation energies _E app are of purely dynamical origin, we 
look for a quantitative relation to PEL structure. A consequence of the rugged nature of MBs is 
that escapes from MBs are not single hops between PEL minima, but complicated multi-minima 
sequences. We introduce the concept of return probabilities to the bottom of MBs in order to judge 
whether the attraction range of a MB was left. We then compute the energy barriers that were 
surmounted. These turn out to be in good agreement with the effective depths i? app , calculated 
from dynamics. Barriers are identified with the help of a new method, which accurately performs a 
descent along the ridge between two minima. A comparison to another method is given. We analyze 
the population of transition regions between minima, called basin borders. No indication for the 
mechanism of diffusion to change around the mode-coupling transition can be found. We discuss 
the question whether the one-dimensional reaction paths connecting two minima are relevant for 
the calculation of reaction rates at the temperatures under study. 



I. INTRODUCTION. 

More than thirty years ago, Goldstein (!]] proposed to 
view a glass-forming system as a point moving in the 
high-dimensional landscape of the potential energy V(x). 
In this framework he suggested to focus onto the local 
minima of the potential energy landscape (PEL) , where 
the system is supposed to be trapped at low enough tem- 
peratures. Via occasional transitions to neighboring min- 
ima the system finally relaxes. Owing to the separation of 
time scales, one is able to describe many features of glass 
formers by properties of only the minima. Stillinger and 
Weber g] formulated this idea in the language of statisti- 
cal thermodynamics using the concept of basins. A basin 
of a given minimum is defined as the set of configurations 
that reach this minimum via their steepest descent path 
x = F{x). (We set x and F(x) as shorthands for the mul- 
tidimensional vectors all particle positions and all forces, 
respectively.) The resulting tiling of configuration space 
into different basins allows one to write the free energy 
approximately as a function of static properties of min- 
ima, i.e. their energies and vibrational frequencies || |J. 
Knowledge of the thermodynamics is in general not suf- 
ficient to predict dynamical properties like diffusion con- 
stants or relaxation times. However, experimental || as 
well as simulated j(| 0, || data seem to indicate that 
there exists a strong connection between dynamics and 
thermodynamics via the Adam-Gibbs relation |1. 

Our goal is to reach a quantitative understanding of 
the slowing down of dynamics, as expressed by the bulk 
long-time diffusion constant D(T). Mode coupling the- 
ory (MCT) jnj predicts a power-law behavior of the form 
1/D(T) oc (T-T c y above the MCT critical temperature 
T c . Since T c is found to be higher than the glass tran- 
sition temperature T g , the MCT divergence of l/D(T) 



at T c is not observed in practice. The common expla- 
nation for this shortcoming of MCT is that the theory 
neglects 'activated processes', or 'hopping', which is sup- 
posed to come into play around and below T c . Indeed 
it was proven by Schr0der et al. jllj that in the vicinity 
of T c the time scale of fast local dynamics around single 
minima becomes well separated from the time scale of in- 
terbasin transitions. Above T c the common picture sug- 
gests that the dynamics are 'entropy-driven' i.e. char- 
acterized by the search for escape directions [|12j [l3| [l4| , 
since saddles lie far below the instantaneous potential 
energy of the system and thus represent no serious bar- 
riers |H| [l6|. Also different observables like the average 
order of saddles |l5, |l7|] and the number of diffusion- 



like normal modes |gj |13|, ^ @ 

seem to indicate that 

well above T c the dynamics is not governed by activated 
transitions between adjacent minima. However, in the 
multi-dimensional space of particle coordinates, it is not 
obvious how to distinguish thermally activated from en- 
tropically limited dynamics. One possible way to do this 
will be discussed below. For the time being, we use the 
term hopping only in a formal sense, meaning that the 
trajectory of the system is mapped onto a sequence of 
jumps among minima. 

Looking for a quantitative link of bulk diffusion D(T) 
to PEL properties, we recently investigated hopping dy- 
namics on the PEL in greater detail |19f . A priori, tem- 
poral and spatial aspects of hopping events have to be 
considered, the former in the shape of the waiting time 
distribution (WTD) of jumps, the latter by the jump 
lengths and directions, and correlations thereof. We 
found that strong backward correlations of jumps arise 
from the organization of minima into superstructures, 
which, following pC| , we call metabasins (MB). It had 
already been known from previous work pl[ that struc- 
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tural relaxation corresponds to jumps among MBs rather 
than single basins. MBs were identified with the help of a 
straightforward algorithm such that close-by minima be- 
tween which the system performs several back- and forth 
jumps arc identified as a single MB. Then, indeed, hop- 
ping among MBs was found to be close to a random walk 
with a distribution of MB waiting times. Motivated by 
this fact, we expressed D(T) in the simple form, 



D(T) 



6N (r(T))' 



(1) 



with the mean waiting time (r(T)) and the effective jump 
length a(T). With this ansatz, we anticipated that wait- 
ing times would carry the major part of the temperature 
dependence. Indeed, a(T) turned out to be constant for 
T < 2T C , which is why we dropped the argument of a(T). 
Eq. [I] constitutes an important step towards the under- 
standing of diffusion in supercooled liquids: it suffices to 
look for the physics behind MB waiting times, spatial 
details of hopping being expressed by a single constant. 

A simple model for hopping dynamics has been dis- 
cussed by Bouchaud and coworkers [p2[ . They consider 
the relaxation from traps of depths E with distribution 
p(E), and escape rates j(E,T) = 70 exp(— fiE). When 
p(E) oc exp(— E/T x ), the WTDs assume power-law tails 
■0(r) cx r" a ( T ) with exponents a(T) = 1 + T/T x . The 
consequence is the divergence of the mean waiting time 
at temperature T x . In our recent paper, we observed 
that the WTDs of a binary Lennard- Jones system are 
in conformance with this kind of power-law decay |L9| |. 
As a consequence of such slowly decaying WTDs, the 
mean value (r(T)) was found to be dominated by the 
few, very long waiting times. In other words, the tem- 
perature dependence of D(T) follows alone from the du- 
rations of trapping in the very stable MBs. These results 
were obtained for small binary Lennard- Jones mixtures of 
N = 65 particles. For a macroscopic system, which, due 
to its dynamic heterogeneity f23|| , contains many slow and 
fast subsystems in parallel, this implies the dominance of 
slow regions in the temperature dependence of D(T). 

The logical continuation along this line of thinking is 
to relate MB lifetimes to the PEL topography. The most 
prominent characteristics of a MB is, of course, its en- 
ergy e MB , which is defined as the lowest energy of all its 
constituent minima. It is then natural to introduce the 
mean MB lifetime (r(e MB ;T)) at constant e MB . Knowl- 
edge of (r(e MB ; T)), together with the population of MBs, 
p(e MB ; T), is sufficient to calculate (r(T)) and thus D(T), 
as we will show now. We write 



(r(T)) 



;T))ip(e MB ;T), 



(2) 



where ip(e MB ; T) is the distribution of MBs visited at tem- 
perature T. We will see that this decomposition can be 
achieved by a detailed analysis of the hopping dynamics. 
Since p(e MB ]T) denotes the probability that at a given 
time the system is in a MB with energy e MB , it is propor- 
tional to ip(e MB ;T) and the time (r(e MB ;T)} the system 



remains in MBs of this energy. With the appropriate 
normalisation one gets 

p(e MB ;T) = — 7— -— <p(e MB ;T). (3) 

From Eqs. [I], || and |^, it immediately follows the repre- 
sentation 



D(T) 



6N \(r(e MB ;T)> 



(4) 



Here, (...) T denotes the canonical time average (w.r.t. 
p(e MB ;T)), while (...) is the average over MBs. Hence, 

{<t( £mb ; T)> , p(e MB ;T)} -> <t(T)> -» D(T), (5) 

where the second implication has been established in our 
recent paper ]is| . The population p(e MB ; T) is related to 
the single-basin population p(e; T), a purely static quan- 
tity, which has been extensively discussed in the litera- 
ture H H fjj. It has turned out for Lennard- Jones mix- 
tures that the number density G eft (e) of minima is ap- 
proximately gaussian. Thus, the population of minima, 
p(e;T) oc G off (e)e _/3e , could be expressed by three pa- 
rameters describing global PEL structure. In the present 
paper, we focus on (r(e MB ;T)), our goal being to deduce 
it from PEL structure. If this succeeds, we have estab- 
lished the following connection, 

local + global PEL structure — > long-time dynamics, 

which, in our opinion, pushes the understanding of diffu- 
sion in supercooled liquids a step further. 

We proceed as follows. We first compute MB lifetimes 
from ordinary simulation, and later compare them to the 
prediction from PEL structure. First, we characterize the 
relaxation from four single, randomly selected MBs. By 
an exhaustive sampling of these MBs, we will be able to 
get some first insights into MB topology. Second, many 
MBs of fixed energy are considered and their lifetimes 
(t^mbST 1 )} are calculated. Third, we relate MB lifetimes 
to PEL structure, by quantifying the MB depths, or ef- 
fective barriers, which determine the temperature depen- 
dence of (r(e MB ;T)}. The physical scenario which will 
emerge from the results of this paper implies that MBs 
can be regarded as traps, surrounded by high barriers. 
It turned out from exhaustive explorations of PEL con- 
nectivity |24| that due to the high dimensionality of con- 
figuration space the number of escape paths from every 
minimum is enormous. Thus, one may anticipate that the 
effective barrier to leave a specific MB results as a com- 
plex superposition of individual escape paths. Therefore, 
enormous numerical effort is required to quantify their 
multitude for many different MBs. 

Note that the whole analysis will be carried out in the 
spirit of activated barrier crossing. The extent to which 
this is present in supercooled liquids is quite disputed 
in literature. However, we will show that for tempera- 
tures in the landscape-influenced regime below 2T C , the 
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apparent activation energy 

E^ p (e MB ;T) = — In (r(e MB ;T)) , (6) 

can indeed be identified with PEL barriers much larger 
than fc^T which the system encounters when leaving a 
MB. Thus, together with Eq. [|, we will find that the acti- 
vated escape out of deep traps is the physical mechanism 
behind diffusion. 

To our knowledge, such a connection between dynam- 
ics and PEL barriers has never been established for a 
fragile glass former. In contrast, for SiC>2, the apparent 
activation energy of diffusion below T c could be related 
to the simple breakage of Si 
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O bonds (25 

The organization of the paper is as follows. In sec- 
tion ||, we provide the details of our simulation, and de- 
scribe the interval bisection method to identify MBs. Sec- 
tion deals with the computation of apparent activa- 
tion energies from relaxation dynamics. The correspond- 
ing energy barriers will be addressed in section after 
introducing our technique for finding transition states 
(section IV). In section VI, we independently demon- 



strate that barriers and associated reaction paths indeed 
govern relaxation. Fina lly, we discussion further as pects 
of our results in section VII and conclude in section VIII . 



finite-size effects ||. Here we choose N = 65, since the 
BMLJ60 system has a stronger tendency to be trapped 
in crystalline configurations. We stress here that the re- 
sults obtained for the BMLJ65 system show no finite-size 
related artifacts. For example, D(T) of the BMLJ1000 
is identical to D(T) of the BMLJ65 above T c , see Fig. |. 
In the temperature range studied, we found that the be- 
havior of a BMLJ130 system largely resembles that of 
two independent copies of a BMLJ65. Thus, the general- 
ization of the present work to larger systems should not 
bear any pitfalls. 

Interestingly, the BMLJ65 relaxation becomes 
Arrhcnius-likc for low temperatures. Since we have 
no equilibrium runs of our BMLJ1000 below T c it is 
unclear whether the Arrhenius behavior of the BMLJ65 
down to T = 0.4 is a finite-size effect. A possible 
explanation is that the lower end of the PEL is reached 
(located at e min w —302, see Fig. |J), preventing further 
increase of barriers. In turn, this may be related to the 
fact that cooperative regions cannot grow any further 
in the BMLJ65, i.e. structural optimization - which 
happens upon cooling - finally comes to an end. Since 
we work above T c , our key results are not affected by 
this argument. 



II. SIMULATION DETAILS. 

A. General. 

In the present work, we investigate a binary mixture 
of Lennard- Jones particles (BMLJ), as recently treated 
by two groups |l7], E?]; see also [p8| . It is characterized 
by the interaction potentials 

V aP {r) = U aP [{a a0 /r) 12 - {a af3 /rf} 

with the parameter set N = Na + Nb = 52 + 13 = 

65, GAB = 0.8CTAA, &BB = 0.88OVL4, CAB = l-S^AA, 

(bb = 0.5e,4yi, r c — 1.8. Linear functions were added 
to the potentials to ensure continuous forces and ener- 
gies at the cutoff r c . These modifications of the original 
potential by Kob and Andersen |2§| are necessary for the 
simulation of small systems. We use Langevin molec- 
ular dynamics simulations (MD) with fixed step size, 
A 2 = 0.015 2 = 2kBTAt/m(, equal particle masses m, 
friction constant £ set to unity, and periodic boundary 
conditions. Units of length, mass, energy, and time are 
oaa, tu, tAAi and m£A /2eaa, respectively. However, 
we will omit these units, for convenience. The mode- 
coupling temperature is T c = 0.45 ± 0.01 in this model 
system (compare p8| ). For the analysis of dynamics from 
the PEL perspective it is essential to use small systems, 
as has been stressed in the literature [^l[ |2{| p0| . On 
the other hand, naturally, the system should not be too 
small in order to avoid major finite-size effects. For the 
BMLJ, N ss 60 turns out to be a very good compro- 
mise JT7|, I2I], p7fl, whereas N < 40 already causes large 



B. Interval bisection. 

By regularly quenching the MD trajectory x(t) to the 
bottom of the basins visited at time t, as proposed by 
Stillinger and Weber, we obtain a discontinous trajec- 
tory A problem from the standpoint of simula- 
tions is to resolve the elementary hopping events. Since 
computer time prohibits to calculate the minimum £(t) 
for every time step t, we normally find ourselves in the 
situation of having equidistant quenched configurations 
U = iAt, with, say, At = 10 5 MD steps. If 
the same minimum is found for times U and tj, we 
need not care about transitions in the meantime, be- 
cause no relaxation has occurred there. If, in contrast, 
7^ we must not expect £(tj+i) to be the di- 
rect successor of £(ij), since many other minima could 
have been visited between £j and fj+i. Therefore, further 
minimizations in this time interval are necessary. For 
reasons of efficiency, we apply a straightforward interval 
bisection method, which locates transitions to an accu- 
racy of 1 MD step: provided £(d°L) + £(da rt ), (a) set 
f(°) <— tf£2rt, *- ^tarti (b) reconstruct the trajectory 
x(t) at time t^ = (t(°) +iW)/2, (c) calculate £(t (2) ), 

(d) if C(t (2) ) = £(£ (0) ), set t(°) «- t( 2 \ else set t« «- 6 2 \ 

(e) repeat (b)-(d) until t^ -t^ = 1 MD step. Repeated 
application of the interval bisection to a simulation run 
x(t) finally gives all relevant transitions. Note that the 
determination of all transitions including the numerous 
recrossings of basin borders would require minimization 
for every MD stepl The interval bisection method thus 
may oversee back- and forth motions between minima 
which, in any event, are irrelevant for relaxation. Al- 
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though computationally demanding, the above method 
has proved most efficient for resolving the relevant de- 
tails of hopping on the PEL and is predestined for the 
construction of metabasins (see below) . 



III. ACTIVATION ENERGIES FROM 
METABASIN LIFETIMES. 

A. Metabasin lifetime construction. 

As said above, stable configurations in the supercooled 
liquid are rarely due to single minima on the PEL, but 
mostly correspond to groups of strongly correlated min- 
ima. While the system is trapped in such a MB for a long 
time, a small number of minima is visited over and over 
again. This is well reflected by the time series of poten- 
tial energies, e(t) — V(£(t)) |h| ^l). In this section, we 
will dwell on the computation of mean MB lifetimes, (i) 
for single, selected MBs, and (ii) averaged over MBs of 
a given energy e MB , thus yielding (7"(e MB ; The indi- 
vidual MBs of (i) correspond to long-lived MBs and thus 
represent typical MBs which govern the temperature de- 
pendence of (t(T)). 

For the grouping of minima into MBs and the resulting 
determination of their lifetimes from a regular simulation 
run, we use the following, straightforward algorithm . 

(a) determine the regions [t*,t\] where t* is the time of 
the first and t\ the time of the last occurrence of mini- 
mum (b) any two regions overlapping by less than 
r mol are cut so that [t*,t\] n = where r mol is a 
small molecular time-scale, (c) any two regions overlap- 
ping by more than r mol are combined to [t|,tj] U [£J,i]i], 
(d) the lifetimes of MBs are defined by the regions after 
step (c), (e) the MB energy e MB is defined as the lowest 
energy of minima visited during the MB lifetime. 

A few comments on the procedure are in order. Time 
regions in (a) are determined by the interval bisection 
method which yields the time of transition from one min- 
imum to another with an accuracy of one MD step. Step 

(b) is motivated by the observation that recrossings of a 
basin border during a transition are very probable. If we 
ignore this fact, i.e. set r mol = 0, step (c) would merge 
nearly all regions and we would end up with unphysically 
long MBs. Instead, we use r mol = 40, which is the basin 
equilibration time obtained from energy autocorrelation. 
Since the durations of transitions are of the same order 
of magnitude as r mol , this is a sensible choice. However, 
we found that the results for MB lifetimes are not very 
susceptible to the precise value of r mol . Step (c) itself 
is the realization of the MB concept. It is important to 
note that, different from plj , we will treat all MBs on 
the same footing here, no matter if they are short-lived 
or long-lived. 

So far, the MB lifetime construction rests upon single 
trajectories, which only partially reflect the configura- 
tion space topology. In section M, the MB concept will 
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FIG. 1: Mean lifetimes of four low-lying, randomly selected 
metabasins, computed from repeated escape runs (e M B = 
-301.64,-300.47,-300.16, and -300.74, from top to bot- 
tom). The number of runs are 85,59,175, and 105, from top 
to bottom. Arrhenius fits work well in the temperature range 
T < 1 ~ 2.2T C , the corresponding activation energies are 
given in the figure. Curves have been shifted vertically by 
0.5(4 — i) orders of magnitude for better inspection. 



be given a more precise, static definition, based on the 
return probability to the ground minimum. 



B. Activation Energies for Single MBs. 

As noted above, the temperature dependence of D(T) 
is dominated by the long-lived MBs. Generally, these are 
low-lying MBs, i.e. deep traps in the PEL. Since different 
MBs differ in their stability, a statistical treatment will 
be needed. As a first step, however, we restrict ourselves 
to the investigation of single MBs. 

The relaxation times computed in this section do not 
stem from regular, linear simulation runs, but are ob- 
tained by artifically placing the system in a specific MB 
and waiting for its escape ('escape runs'). The above 
algorithm for the MB lifetime construction implicitly as- 
sumes that MBs finally have been left. In other words, 
the algorithm may not be used to determine the time 
where to stop the simulation due to successful escape. 
Fortunately, we can avoid running into this paradoxi- 
cal situation by judging from an independent criterion 
whether an escape has been completed: if the distance 
of the instantaneous minimum to the starting position is 
greater than d max = 4, returning to the original basin can 
practically be excluded (see section |y| for a justification 
of <i max = 4). Then, by applying the MB construction al- 
gorithm to the escape run, we obtain the lifetime of the 
MB. 

We analyzed four low-lying (e M B < —300), randomly 
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selected MBs in greater detail. By repeated starts from 
the bottom of the MBs, we computed the mean lifetimes 
(t,-(T)) as a function of temperature. From Fig. [I], we 
see that the relaxations from all MBs follow nicely an 
Arrhenius law. We note that, due to starting at minima, 
a short intra-basin equilibration time (r mol = 40, from 
energy autocorrelation) has been subtracted from the raw 

(n(T)). 

The fact that an Arrhenius form of (T(e MB ;T)) is ob- 
served indicates that the barriers do not change any fur- 
ther upon lowering temperature. Put differently, MBs 
serve as traps surrounded by barriers with heights around 
E app (i) = din (r;(T)) /d/3. We will see in section that 
this is indeed correct. Since E app (i)/kBT c > 10, this im- 
plies a strongly activated dynamics near T c . 



C. Activation Energies vs. MB energies. 

As a further step, we analyze the mean relaxation time 
from MBs with the same energy, (r(e MB ;T)); see Eq. [|. 
Clearly, the low e MB are not populated at high tempera- 
tures so that regular simulation does not yield (r(e M B; T)) 
over a wide temperature range. We therefore artificially 
place the system in the desired MBs (in the lowest min- 
ima e M B thereof) and measure the escape times as a func- 
tion of temperature. Averaging over many different MBs, 
we obtain (r(e MB ;T)). Results are shown in Fig. || as 
a function of e MB . Below T = 1, all relaxation times 
display Arrhenius behavior. Thus the apparent activa- 
tion energies i? app (e MB ;T) are temperature independent. 
In the following we will therefore omit the second argu- 
ment. Thus, we can write 



(r(e MB ;T)) = r (e MB )e^^\ 



(7) 



As expected, the properties of MBs as expressed by 
B app (e MB ) depend on their ground state energy e MB . 

We can interpret -f? app (e MB ) as the mean effective depth 
of MBs at e MB . Since the lower end of the energy land- 
scape is reached at e « —302 no deeper traps exist (com- 
pare Fig. ||, see also (50j). A simple statement for the 
depths of traps would follow if the rims of all traps 
were at the same level e th . The consequence would be 
E &pp (e UB ) = e th - e MB , for all e MB < e th . This simple sce- 
nario is ruled out by the data, see Fig. ||(b). Actually, 
a more complicated energy dependence of -E app (e M B) is 
expected from the very fact that the system - despite its 
small size - is not a comp letely cooperative unit, see the 
discussion in section VI] . 



The fact that we still observe Arrhenius-like relaxation 
in Fig. U indicates that the variation of trap depths at 
constant e MB is not large, compare E app (i) from Fig. [j]. 
Otherwise, E app (e MB ) would increase upon decreasing 
temperature, due to the more and more dominant, ex- 
tremely deep traps. In contrast, trap depths at constant 
e MB seem rather well defined by e MB , which suggests the 
existence of some underlying topological principle. 




0.0 0.5 1.0 1.5 2.0 2.5 
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log -C (£mb) 
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-300 -295 -290 -285 -280 -300 -295 -290 -285 -280 

FIG. 2: (a) Arrhenius plot of mean MB lifetimes (r(e M B;r)), 
for different £mb • A MB equilibration time of r mo i = 40 was 
subtracted. Straight lines are fits of the form Eq. [?]. (b) ap- 
parent activation energies £a PP (eMB). (c) prefactors to(£ M b). 
Curved lines are interpolations of the data. 



As seen from Fig. ||(c), the prefactor r (e M B) has no 
strong dependence on Cmb- From high energies, it de- 
creases at most an order of magnitude and seems to level 
off below e M B = —297. Hence, for the range of ener- 
gies that dominate (r(T)) at low temperatures, it can 
be considered constant within error bars. In contrast to 
E B , pp (e MB ), we will not be able to deduce To(e MB ) from 
PEL structure. Its weak variation is therefore quite for- 
tunate. 

We will now analyze the second factor of the integrand 
in Eq. ||, ip(e MB -.T). It is shown in Fig. ||[ Interestingly, 
the variation of ip(e M B',T) is much weaker for low T as 
the variation of p(e MB ;T). From Eqs. || and ^, one con- 
cludes that the constancy of the distribution <p(e MB ;T) 
is equivalent to having -E app (e M B) = e«, — £mb, with some 
constant e th . Since this simple behavior is not present, 
one must still have a residual temperature dependence of 

¥>(£mb; t). 

Concerning p(e MB ; T) , the weak temperature depen- 
dence of its first moment for the three lowest tempera- 
tures is simply related to the probing of the lower end of 
the PEL. Actually it turns out that p{e mB ;T) is, within 
statistical error, identical to the corresponding distribu- 
tion of minima p(e;T). One would expect this for high 
e MB , because no pronounced MBs are observed there. 
Considering a deep MB with many minima, this will 
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FIG. 4: Arrhenius plot of the mean waiting time (t(T)} versus 
the indirectly determined counterpart, (r(T)) i .. For compar- 
ison, we also show the inverse one-particle diffusion constant 
1/D(T) multiplied by a constant (a 2 = 1.0), see Q. Error 
bars are of the order of the symbol size. Also included is the 
1/D(T) for the BMLJ1000 system. 
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FIG. 3: (a) distribution ip(e M B',T) of MB energies, for four 
temperatures, (b) Mean energies, from </j(6mb;T) and from 
p(eiAB',T). (c) Variances of the distributions cp and p. Poly- 
nomial fits to the data are shown in (b) and (c). Straight 
lines are predictions for p from an ideally gaussian distribu- 
tion G e ff(e M B) (mean eo, variance a 2 ). Deviations are due to 
reaching the lower end of the PEL, i.e. the deepest amorphous 
minima. 



equally effect no large difference between p(e MB ; T) and 
p(e;T). The reason is that the group of minima near 
e MB carry the largest part of the population. Since they 
are close to e MB , transferring their weight to e MB when 
computing p(e MB ;T) has little effect. 

As a consistency check, we use the data from Fig. || 
and H to reproduce (r(T)) indirectly via Eq. || (denoted 
(t(Tj). ). The match with (t(T)) is not completely triv- 
ial since the data for (r(T)) and ip(e MB ;T) were gathered 
from a linear simulation run, while (r(e MB ;T)) results 
from selected MBs of certain e MB , where the system has 
been artificially placed. As shown in Fig. [|, the agree- 
ment of (t(T)) and (r(T)) lnd is good for T < 1 within the 
possible accuracy. Note that there is no free fit parameter 
between them. The deviation at T = 2 can be explained 
by the fact that (r(e MB ; T)), above T = 1, and especially 
for the high e MB , departs from Arrhenius behavior, see 
Fig. |(a). 

So far, all barriers or trap depths have been derived 
indirectly, from the temperature dependence of waiting 
times. A link to the PEL structure is still lacking. For 
instance, the activation energies _E app (i) of this section are 
expected to reflect the local topography of the selected 



MBs. Indeed, they can be identified from the barriers of 
escape paths, as will be demonstrated in section 

First of all, the barriers between neighboring minima 
are of great interest. These arc known once we have in 
hand the corresponding transition states. 



IV. NON-LOCAL RIDGE METHOD FOR 
FINDING TRANSITION STATES. 

A. Description of the method. 

We now describe how to determine transition states 
(TS) from the simulation, by what we call the (non- 
local) ridge method. The principle idea is that TSs are 
local minima of basin borders. They can be pictured 
as the lowest points of mountain ridges on the PEL. If 
the system crosses a basin border at time t, the steep- 
est descent path starting from x(t) should end up in a 
TS, see In practice, however, the descent will de- 

viate from the ridge due to numerical error, finally end- 
ing up in the minimum £o = £,(t— ) or £i = As 
a way out, we let the system perform two descents in 
parallel, on either side of the basin border, as schemati- 
cally depicted in Fig. ||. More specifically, if a transition 
happened at time t, interval bisection yields the config- 
urations x(t) = y and x(t + 1 MD step) = y\. From 
these, by further interval bisection on the straight line 
between yo and yi, the distance to the border may be 
further reduced if necessary, resulting in two configura- 
tions, again called yo and y±. Close as they are, they still 
belong to different basins. If we now let descend yo and 
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FIG. 5: Sketch of the TS search with the ridge method. 



yi in parallel, they first move along the ridge towards 
the transition state until they finally bend off to their re- 
spective minima. This separation is clearly not wanted, 
so from time to time we reduce their distance by inter- 
val bisection. After a few iterations (descents+interval 
bisection) the vicinity of the transition state is reached 
in most cases. We then use a short minimization of the 
auxiliary potential V = \ \F(x)\ 2 followed by a few steps 
of Newton-Raphson type, which bring the search for the 
TS to a quick convergence. Besides a vanishing force, the 
resulting configuration £ has a Hessian matrix with one 
negative eigenvalue. After small displacements along the 
corresponding eigenvector, one reaches the adjacent min- 
ima via steepest descent. This yields the reaction path 
(RP) C(s), where s is a curvilinear parameter. We set 
C(0) = Ci C( s o) = Co, and C( s i) — £ii where either sq or 
si is negative. 

It can happen, though, that no saddle between yo and 
yi is found, but that the interval bisection locates a third 
minimum. The basin border splits into two at this point, 
and no direct saddle between the initial and final mini- 
mum is available. Thus, we also have to split the descent 
along the basin border into two processes and then con- 
tinue separately If the two descents are successful with- 
out further bifurcations, we are finished and have the 
optimum reaction path which takes a detour via a third 
minimum. In such a situation, the RP is clearly not very 
useful. It has to be stressed that bifurcations are no ar- 
tifacts of the ridge method, but a topological feature of 
some basin borders on the PEL. Fortunately, as a sig- 
nature of strong anharmonicity, they are quite rare and 
happen to occur only in the high-energetic regions of the 
PEL. For the escapes from long-lived MBs, they are of 
no importance. 

A similar algorithm is described in the literature 
which, instead of minimization and interval bisection, 
uses local maximization between yo and y\ to prevent 
the configurations from moving apart. Although compu- 



tationally less expensive, this method is not appropriate 
for our purpose. As an effect of the high dimensional- 
ity, the local shape of the PEL around yo and y\ gives 
no direct clue to the membership to basins. When de- 
scending, one may thus loose the important property of 
yo belonging to the basin of £o and yi belonging to that 
of £1 . This effect has indeed been reported in |n] . 

In the literature, plenty of methods exist dealing with 
the computation of transition states. One kind of them 
starts from the knowledge of the initial and final mini- 
mum (|^, H||, |35| . After a more or less educated guess 
for an initial trial RP, one iteratively improves the RP 
according to some prescription, e.g., the minimization of 
an action functional. Two sources of erroneous results 
have to be addressed in this connection. First, the two 
minima in question have to be true neighbours. This can 
only be verified by locating two points close to the basin 
border, e.g. by interval bisection of the initial trial path. 
The numerical cost is not small; for our ridge method, for 
instance, about one third of the calculation time is con- 
sumed by fixing yo and y\ (depending on the minimiza- 
tion interval of the MD run). Second, the iterative path 
optimization may become stuck in a local extremum, due 
to an unfortunate choice of the initial path. 

The other kind of TS search methods start from an ini- 
tial minimum and climb up to a transition state guided 
by the shape of the PEL. Just walking against the force, 
however, would be a fatal strategy, as one can see by 
turning the PEL upside down: ending up in a TS is 
numerically impossible, since one quickly runs into one 
of the PEL singularities (two or more identical particle 
positions). Eigenvector-following algorithms |36| over- 
come this defocussing of steepest ascent paths by walk- 
ing into the direction of negative local PEL curvature. 
The 'activation-relaxation technique' by Mousseau and 
coworkers, in contrast, steps against the force in the 
direction leading away from the minimum, while de- 
scending the PEL perpendicular to that direction [p7| . 
A drawback of the latter methods is that the choice 
for the next TS to mount is not well under control. 
From the minimum, a starting direction is chosen, ei- 
ther by purely random displacements or by some hard- 
sphere- like particle moves pi] . Unfortunately, the num- 
ber of escape directions from a minimum is generally 
very large (at least O(Nd) as we found in the BMLJ65), 
whereas the majority of these is dynamically inaccessible 
at low T. Hence, eigenvector-following and activation- 
relaxation techniques yield many TSs which only negli- 
gibly contribute to relaxation rates. Striving for the sim- 
ulation of low-tem per ature hopping dynamics based on 
these methods |^7|, |9| , one may suffer a considerable 
reduction of efficiency. 

We finally mention two complementary means of 
studying energy barriers. The 'lid' algorithm, proposed 
by Schon and coworkers [ [l0| , is able to find upper bounds 
for the depths of single basins. By performing random 
walks below different potential energy thresholds and by 
regular minimizations, one is able to compute the ele- 
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FIG. 6: Comparison of transition states, obtained via the 
ridge method, with minima of the auxiliary potential V. 
Starting points for saddle computations lay close to basin bor- 
ders. Main plot: histogram of t^-saddle minus TS energies. 
Inset: histogram of indices of V-saddles. 



vation necessary for transitions to neighboring minima. 
From a more theoretical perspective, Schulz has speci- 
fied a relation between transition rates and the overlap 
of vibrations in neighbouring basins EL 



B. Comparison to V-saddles. 

The advantage of the ridge method is that we defi- 
nitely find the relevant barrier for a transition, i.e. a 
first order saddle on the basin border next to the point 
where the border was crossed. In contrast, the method 
using the auxiliary potential V = ^\F(x)\ 2 as applied 
in recent studies jl5|, [16], [l7|, ^5(| has two major draw- 
backs. First, the V minimization locates saddles (we 
call them I^-saddles), even if they are not accessible ki- 
netically. This is because the expression F^HF is not 
positive (H = H(x) denotes the Hessian of V(x)), i.e. 
^-minimization can climb up to a saddle. Second, one 
obtains higher-order saddles and, most frequently, non- 
stationary points (shoulders) . These configurations are of 
no use to us because we specifically analyze paths over the 
lowest barriers on basin borders, i.e., transition states. 

To shed more light on the interrelation of TSs and V- 
minima, we minimized V by steepest descent, starting 
from configurations x(t) only if £(t) ^ £(t + 1 MD step) 
(like yo in Fig. ||). In other words, we calculated V- 
saddles exactly at transition times. If this yielded the cor- 
rect TSs, our more time-consuming ridge method would 
be clearly useless. The difference Ae = e v — e T s specifies 
the overestimation of the true barrier by the V^-saddle. It 
may also happen that the index of the I^-saddle is differ- 
ent from one. The distributions of Ae and the index are 
shown in Fig. ^| (T = 0.5). Obviously, the V^-saddles con- 
siderably overestimate barriers and the correct TSs are 
only found very rarely. Moreover, most of the V^-saddles 
have an index different from one, i.e. are no TSs at all. 



In turn, the energy of the TS is never undersold by a 
F-saddle. In conclusion, V^-saddles turn out to have the 
undesired quality of being decorrelated from the relevant 
TSs, i.e., from the barriers that control relaxation (see 
section [v]) . 



C. Population of Basin Borders. 

After Angelani and coworkers Jl5|, [l(J, the mean in- 
dex of ^-saddles vanishes at T c . Therefore, as they have 
argued, dynamics above T c is dominated by saddles, in 
that there are always some unstable directions available 
which allow the system to relax, without traversing an 
additional energy barrier. Passing T c , the mechanism 
suffers a drastic change, and abruptly, one is faced with 
an index of ca. zero, i.e., saddles have to be reached 
via thermal activation. Since the preceeding subsection 
may cast some doubts on the significance of ^-saddles, 
we now want to discuss an alternative analysis of the way 
the population of minima versus unstable configurations 
evolves upon decreasing temperature. More specifically, 
we determine the population of basin borders, 



Z(T) 



dB / dxe-^^Six-B), (8) 



where integration is over the non-crystalline part of con- 
figuration space, also in the partition function Z(T), and 
B runs over all basin borders of the PEL. This expression 
is impractical in numerical simulation; one may rather 
ask if, for some instantaneous configuration x, there is 
a basin border nearby. In this case, small random dis- 



placements (length 8 £ IR, direction lj G ]R 



JY,/ 



1) 



possibly lead into another basin, i.e. ^ £(x + u>8). 

This kind of PEL analysis has been recently carried out 
by Fabricius and Stariolo J42|. One calculates 



Pbb (TiS) = (P(Z(x)?S(x + w8))) 



T.ui ' 



(9) 



which is the probability that random disturbances u>8 will 
cause crossings of basin borders at temperature T. The 
brackets denote the canonical plus the average over the 
random directions uj. One obtains the behavior 



p BB (T; 8) -> const x p BB (T)S, 8^0 



(10) 



(the constant is set to unity for convenience). The valid- 
ity of Eq. [l0| is demonstrated in the left inset of Fig. 
where p BB (T; 8)/ 8 has been calculated as a function of 
We find that p BB (T; 8)/ 8 is constant within statistical er- 
ror below 8 = 1.2. As an orientation, the typical distance 
between neighboring minima is larger than 2.0, whereas 
intra-MB neighbors on average are less than 1.0 apart. 

The main part of Fig. shows results for p BB (T) in 
an Arrhcnius plot, with 8 = 0.7. Over the whole tem- 
perature range considered, p BB (T) is Arrhenius-like. The 
apparent activation energy is ca. 1.8, which is small in 
comparison with the typical values observed for MB life- 
times. However, the temperature dependence becomes 
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FIG. 7: Population of basin borders, p BB (T), obtained from 
disturbances of length S — 0.7, which corresponds to a dis- 
placement of ca. 0.09 per particle. Left inset: dependence on 
5 of Pbb(T), for T = 0.5 and T = 0.6. Right inset: p BB (T) 
plotted linearly against T. 
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FIG. 8: Potential energy along the reaction path which 
was calculated from the dynamics during 10 5 MD Steps, at 
the end of a typical MB of life span 8 x 10 6 MD Steps. The 
mapping of s to time is non-linear. The small barriers for 
s < 5 belong to fast intra-MB transitions. p ba ck denotes the 
probability of returning to the bottom of the MB. As a com- 
parison, the potential energy at that temperature (T = 0.5) 
fluctuates around —249.3 ± 6.1. 



stronger if we impose the constraint of a minimum dis- 
tance between neighbouring minima (data not shown). 
In this way, we eliminate the fast intra-MB transitions, 
which have small barriers. 

In any event, p BB (T) features no noticeable change in 
behavior when approaching and crossing T c . In a differ- 
ent graphical representation (see right inset) one might 
wrongly conclude that p BB (T) disappears at some finite 
temperature. Stated differently, the data suggest that 
the increasing timescale separation upon cooling happens 
rather smoothly, with no distinctly new physics emerging 
near T c . This is in qualitative agreement with the work 
of Schr0der et al. [ri|, who use the incoherent scattering 
functions from hopping dynamics £(i) to deal with the 
separation of intra- and interbasin dynamics. There, the 
initial short-time decay of scattering functions (quanti- 
fied by the so-called non-ergodicity parameter) is nothing 
else than a measure for the population of basin borders. 



V. ENERGY BARRIERS FROM PEL 
TOPOLOGY. 

A. Return Probabilities and Metabasin Definition. 

With the tools of interval bisection and TS search, we 
are now in the position to analyze the escapes from MBs 
in full detail. When a MB is left, we first resolve all min- 
ima visited during the escape. Second, all corresponding 
TSs and, if desired, reaction paths are calculated. An ex- 
ample is shown in Fig. ||. The successive RPs were spliced 
together to a long, multi- minima RP C( s )- One might 
take the energy profile, VX£(s)), depicted in the figure, 
for one of the common cartoons of a PEL. However, it 



rests upon real data. Berry and coworkers have produced 
similar charts for the relaxation of small atomic clusters 
towards their global minima |39], For s < 5 one can 
sec the typical back- and- forth hopping among the ground 
minima of the MB. Obviously, the corresponding barriers 
are not large compared to k-gT = 0.5. The escape starts 
at s = 5. The first minimum reached is very unstable as 
expected from the small backward barrier. Indeed, if we 
repeatedly start in this minimum and perform a number 
of short simulation runs (here: 99) with different random 
numbers, the system will return to the bottom of the MB 
with probability p back — 98% and leave the range of at- 
traction only rarely. Thus, the escape is far from being 
complete at this stage. Going to the next minimum, the 
return probability decreases, but does not drop to zero. 
We say that the system is free if p back is smaller than 
50%. As the outcome of this investigation, we obtain 
the energy barrier surmounted before the minimum with 
Pback < 50% was reached, see below. The exits from other 
long-lived MBs mostly look the same as in the example, 
while the escape in one jump is not common. In other 
words, MBs usually have the form of a funnel with some 
ledges on the walls |2(| Minima with p back > 50% 
are said to belong to the MB. This criterion is reminis- 
cent of the definition of dynamic bottlenecks introduced 
by Chandler and coworkers |Q . 

An interesting property of a MB is its diameter d. It 
is defined as the maximum distance between its minima. 
For the MBs found in the simulation at T = 0.5, the 
distribution of diameters is depicted in Fig. ||. The delta- 
peak from single-minimum MBs has been omitted. No 
MB with d > d max = 4 has been found. As a consequence, 
if a minimum has a distance larger than ri max to some 
MB minimum, we can safely assume p b ack *C 50%. This 
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FIG. 9: Distribution of MB diameters, d, defined as the max- 
imum distance between all minima that were visited during a 
MB lifetime. The delta-peak from single-minimum MBs has 
been omitted. 




Kernel 



Pback>50% 



FIG. 10: Sketch of the configuration space around a MB, 
crosses representing minima. Large crosses are the highly 
populated minima on the bottom of the MB. The shaded 
area comprises minima of high return probability to the kernel 
minimum (p back > 50%). By definition, these constitute the 
MB. The bent line is the system trajectory x(t) entering and 
finally leaving the MB. 



on temperature, since it is defined by dynamics. Corre- 
lations among minima are expected to increase towards 
lower temperatures, implying that MBs are no static 
concept but rather grow with decreasing T. In Fig. ||, 
e.g., the minimum at s ~ 6.5 has the 'critical' value of 
Pback ~ 47% at T = 0.5. Although we do not know the de- 
tails of PEL connectivity around this minimum, the small 
backward barrier suggests that the minimum would ex- 
ceed p back = 50% for still lower temperatures, thus joining 
the MB. However, we may also conceive some situations 
where a critical p back « 50% is quite unsusceptible to 
temperature changes. This is the case if backward and 
forward barriers are of about the same size. We will come 
back to that issue later. 

We further note that the explicit computation of p back 
can be extremely expensive. This is mainly the case when 
p back is small, and complete escapes beyond d max have 
to be awaited. However, the exact value of p back is of 
no great interest. In fact, it suffices to know whether 
Pwk < 50% or p back > 50%. This decision can often be 
reached to a high confidence with few trials. 

The MB lifetime algorithm in section III is based on 



the detection of back-and-forth jumps between minima. 
One mostly observes the dominant minima on the bottom 
of the MBs, whereas the more elevated members are only 
weakly populated, see Fig. [l(]. If MB lifetimes are to be 
read from a simulation run, it suffices to notice when the 
set of dominant MB minima has been left, since the visits 
to the elevated minima at the end of the MB lifetime 
happen quite rapidly. Thus, the algorithm of section III 
reduces the MB to the most populated minima, which is 
sufficient for the purpose of lifetime calculation from a 
given simulation run. In contrast, for the prediction of 
MB relaxation behavior as pursued in this section, the 
minima close to the rim of MBs arc of special interest. 
Their elevations from the bottom of the MB give the 
depth of the MB. 



B. Barriers for Metabasin Relaxations. 



criterion has already been used in section III 



Based on these insights, we can now provide a more 
complete description of MBs (Fig. |l(J). First, the ground 
state of a MB has to be identified (kernel minimum), 
since the definition of p back rests upon it. At low enough 
temperatures, the kernel minimum will certainly be vis- 
ited during the MB lifetime, due to the very low barriers 
among the minima on the bottom of the MB. Second, 
for minima beyond the distance <i max from the kernel, we 
set p back to zero. Third, the probability p back for return- 
ing to the kernel before reaching a distance greater than 
d max can be assigned to every remaining minimum and 
-in principle- be computed by simulation. To this end, 
one repeatedly starts in the minimum and checks if a re- 
currence to the kernel occurs. Fourth, the minima with 
Pwk > 50% are defined as the MB. 

Please bear in mind that p back will in general depend 



In the spirit of the above remarks, we will now carry 
out a systematic investigation of the energy barriers over- 
come when escaping MBs. The goal to is to recover the 
apparent activation energies computed in section III from 
PEL topology. 

The mean lifetime (rj) of MB i can be expressed in 
terms of escape rates 7j iQ of different relaxation chan- 
nels a, 

a 

In general, each 7j iQ reflects a multi-minima escape path 



> Coi f. 1,12 f 

4o — > si — ► 4: 



Ci 



£m-: 



C.M-1 



M 



(12) 



as the one shown in Fig. |g. Here, £o is the kernel min- 
imum and C Q b is the TS for £ a — > Suppose that 
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the number M of jumps in the sequence Eq. [L2| is large 
enough to completely quit the MB's range of attraction, 
i.e., p back (M) ~ 0. For the escape shown in Fig. |L e.g., 
M > 7 would be fine. 

We further take for granted that the rates for sin- 
gle barrier crossings follow quantitatively -via transi- 
tion state theory- from the height of barriers, E a b = 
V((ab) — V(£ a ) (the energy difference between the mini- 
mum £ a and the TS between a and b). Hence, rates g a b 
for single transitions £ Q — > £b are characterized by 



g a b oc e 



-/3E a 



(13) 



A justification for this assumption, even for temperatures 



above T c , will be given in section VI 



Generally, the probability of upward jumps is small at 
low T. Hence, climbing out of a MB in a back-and-forth 
fashion (e.g., £ a = £ a+2 and £ a+1 = £ Q+3 ) is not probable. 
(This is reminiscent of the fact that the activated crossing 
of single potential barriers happens on a short time scale, 
i.e. in a rather straight way.) In contrast, excursions 
from the main path may happen. As shown in Fig. [|, 
the minimum at s = 6.5 is revisited at s = 9 after taking 
a look at another minimum (s w 8). The latter does 
not appear again later on. Clearly, running into such 
'dead ends' should not contribute to the escape rate via 
the successful main path. We therefore eliminate such 
excursions from the sequence of minima, Eq. [l^. From 
these remarks we take the liberty of assuming that no 
minimum appears more than once along the escape path, 



One possibility for calculating activation energies from 
Eq. [l(] would be to consider the complete paths Eq. [l^ , 
where p b k (Af ) « 0, and determine all terms in the sum 
of Eq. |l6|. However, an accurate computation of all the 
desired p re t(a)'s would even be more costly than the de- 
termination of the point where p back changes from above 
to below 50%. We therefore use the following approxi- 
mation of Eq. [l^, which is in conformance with our pre- 
vious definition of MBs: Let m(T) be the first minimum 
along the path Eq. Q where p back < 50%. Then, for all 
a < m(T), we set p lct (a) to unity, while for a > m(T) 
(i.e. outside the MB), we let p ret (a) = 0. Thus, 



_d_ 

~d/3 



m— 1 

ln7^ Q Ei }a = Eq\ + (E a ,a+1 — E a ,a-l) 
a=l 

^m — 1 ^0 ~t~ Ejjyi rn—Xi 



(17) 

where m — m(T). In this way, the terms a < m(T) in 
Eq. HI are given higher weights, whereas those of a > 
m(T) are neglected. We will dwell on the quality of this 
approximation later on. 

Note that, due to the temperature dependence of p back , 
energy barriers Ei >a generally increase upon cooling: At 
high temperatures, in contrast, correlations among min- 
ima are small, such that MBs (even the low-lying) consist 
of only one minimum. This effect is included in Eq. ^ 
by the temperature dependence of m(T). 



£,a 7^6, a ^ b. 



(14) 



We are now interested in the contribution of the path 
Eq. [l2] to the total escape rate Eq. 111]. Particularly, we 
have to consider the question of how many single transi- 
tions are relevant for the escape process. The probability 
to jump from minimum £ a to £ a+i is g a ,a+i/g a , where g a 
denotes the inverse lifetime of minimum £ a . The rate of 
escape via a longer pathway now is given by the rate of 
the first jump times the probability that the minima £ a 
(o = l, ...M) are visited in correct order thereafter, 
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With the help of Eq. O one calculates 



(15) 



M-l 



E01 + > _ p^(a)(E a ^ a+1 - £ 0j0 _i), (16) 



a=l 



where p Iot (a) — g a ,a~i/9a is the probability to jump back 
to minimum a — 1 from minimum a. In the deriva- 
tion of Eq. (l6|, we have neglected a term proportional to 
Ea^a+x minus the average barrier when jumping from a 
to a neighbouring minimum other than a — 1. This term 
strictly vanishes when performing the final summation in 
Eq. |ll|. Moreover, we made use of Eq. [[§ 



C. Single Metabasins. 

We n ow re late the lifetimes of single, selected MBs (cf. 
section [II B) to PEL barriers. By repeated starts from 
these MBs, the local PEL topography is sampled thor- 
oughly, yielding sets of typical escape pathways. When- 
ever a MB is left, we locate the transitions by interval bi- 
section and obtain the corresponding TSs with the help 
of the ridge method. Then, p back is calculated for the 
minima visited, until for the first time, p back < 50%. Fi- 
nally, the barrier E i a ^ is computed according to Eq. [l7| , 
where a(k) denotes the escape path chosen at the fcth es- 
cape. The histograms of barriers are shown in Fig. 
for the four MBs of Fig. at T = 0.5 = 1.1T C . Due to 
the slow dynamics at this temperature, the computation 
of p back was rather expensive. Nevertheless, the statistics 
should be sufficient for a reasonable estimate of the ap- 
parent activation energy. To this end, we express E app (i) 
of MB i in terms of the contributions Ei a , 
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Pi,a.Ei a 



£^,(i), 



_ _ (18) 
where Eqs. [Ll] and [17| were used. Thus, the barriers Ei_ a 

are weighted by the probabilities Pi >a = 7i, Q / X) Q 7i,a 
that the escape happens via pathway a. Note that the 
Ei,a(k) correspond to the pathways that were chosen by 
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FIG. 11: Bold curves: Histograms of barriers E i>a (k) over- 
come when escaping single MBs (i = 1,2,3,4 at T — 0.5). 
Light curves: Respective histograms of barriers Eoi from first 
jumps. Apparent activation energies i? app (i), mean barriers 
i?° PP (i), and mean barriers from first jumps Eoi are given in 
the figure. 



the system, i.e. they are already weighted correctly by 
Pi,a(k), compare Eq. [l8| Therefore, E°^ p (i) is just the 
average of the The values of i? app («) and E^ v {i), 

given in Fig. [ll], are in good agreement. Also shown in 
Fig. [ll] is the distribution of first barriers, Eqi, belong- 
ing to the step £o — > £i. Evidently, the neglect of the 
multi-minima nature of escapes leads to a considerable 
underestimation of apparent activation energies. 

We now continue the discussion of the temperature de- 
pendence of barriers Ei^ a (T). At the example of MB 1 
from Fig. [ll], we have carried out the above program 
for two other temperatures, T = 0.6 and 0.8. The ob- 
tained distributions of barriers, P(£'i jQ ,), are shown in 
Fig. |l2|. We find that the estimates for the apparent 
activation energy (£^(1) = 6.9 ± 0.5, T = 0.6, and 
E^ = 6.8±0.5, T = 0.8) remain in good agreement with 
-S ap p(l) = 6.7 ± 0.3 from section III. The distributions 
of barriers, however, grow narrower with decreasing tem- 
perature. Single, high barriers, contributing to the right 
wing of the distribution, become inaccessible at low T, 
i.e., the relative weights p^ a of the corresponding escapes 
become small. This suppression of high barriers at low 
T is a trivial effect. 

More interesting is the vanishing of small barriers upon 
cooling, i.e., of the barriers E < 5 in the figure. Naively, 
one would expect these to dominate the escape rate at low 
T. However, due to the stronger backward correlations 
(increased p ba ck), jumps over these barriers eventually do 
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FIG. 12: Normalized histograms of barriers E\ } u overcome 
when escaping MB 1, for T — 0.5,0.6, and 0.8. The number 
of contributing barriers are 42, 72, and 59, respectively. Es- 
timated apparent activation energies, E^ p (i;T), are given in 
the figure. 
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FIG. 13: Histogram of barriers from a regular MD run at T = 
0.5 (bold). Neglecting the contributions of the last transition 
state (e m — £ M b), we find smaller barriers (light line). The 
barriers Eoi from only the first jumps are given as the dotted 
line. 



not suffice anymore to escape. As described above, the 
respective escape paths, £o - ► ■■■ - * ?m(T)j grow longer, 
and the barriers change to a different, mostly larger value. 



D. Average over Metabasins. 



During our analysis of the escape times in section III 
the apparent activation energies -E app (e M B) emerged as 
useful quantities. Although the above results already 
indicate that barrier hopping is the relevant motional 
mechanism, a clear-cut verification requires the compari- 
son with the average barrier the system has to cross when 
leaving a MB with energy e MB . 
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For this purpose we now carry out a similar program 
as before, with many MBs visited during an ordinary 
MD run. We concentrate on MBs with lifetimes of more 
than 10 5 MD steps (179 MBs) at T = 0.5. When such 
a MB is left, we locate the transitions by interval bi- 
section and obtain the corresponding TSs by the ridge 
method. Then, we calculate p b ack and identify the barrier 
Ek = Ei(k),a(k) according to Eq. [jj]. The histogram of 
barriers is shown as the bold line in Fig. [l^. For compar- 
ison, we also show the barriers minus the contribution of 
the TSs Er m _x\ m . Ignoring multi-minima correlations, 
we further show the histogram of first barriers Eqi of 
escapes. Evidently, the neglect of TSs or of backward 
correlations leads to much smaller barriers. 

From the above barriers, we will now calculate esti- 
mates of apparent activation energies. When the average 
over lifetimes of different MBs is considered, each MB i 
acquires a weight (fi corresponding to its probability of 
occurrence, 



pest / \ 



Eapp( £ MB) 




Eapp( £ MB) 



(t) =^2 Pi in) 



slope=-0.55 

290 



FIG. 14: E a pp (emb ) (Fig. |) vs. estimated ^(cmb) from 
PEL barriers. Considering only the first jumps of escapes, we 
find a much smaller estimate (^(emb)). Data stem from 
a regular MD run at T = 0.5, where MBs of lifetime greater 
than 10 5 MD steps were analyzed (179 MBs, see Fig. fl3|). 



At fixed e M B , the analog to Eq. [18] can then be derived 

where summation goes over MBs of energy e M B- As in 
Eq. [y| the barriers in Eq. [n] are weighted according to 
their probability of occurrence, but, additionally, with 
the respective MB lifetimes. 

In Eq. [l9|, we have neglected terms stemming from 
the variation of tpi's with temperature. This is justified, 
since the ipi's belong to the same e M B- Their relative 
weights will only vary if these MBs differ considerably in 
barrier heights. As already stated above, however, MBs 
of the same energy seem to be fairly uniform regarding 
this property. For the finite sample of MBs visited during 
an MD run, Eq. [l^ then takes the form 



T kEk 



(20) 



where summation goes over MBs of energy e M B- Again, 
the correct weighting is implicit here. This expression can 
be shown to converge to the right-hand side of Eq. [l9| in 
the limit of infinitely long sampling. In 
the values of E^ p (e un ), determined in this way. 



14 we show 
They 

perfectly agree with the apparent activation energies, de- 
rived from the analysis of relaxation times at different 
temperatures. Thus we have a clear-cut proof that the 
apparent activation energies E app (e MB ) are indeed related 
to barriers on the PEL and thus reflect activated behav- 
ior significantly above T c . This again demonstrates that 
we not only deal with the right order of barrier sizes, but 
we also quantitatively link PEL topography to dynamics. 

For comparison, we included the apparent activation 
energy which results, if only the first transitions of es- 
capes, £o — ► £i, are considered (E i — V^(Coi) — £o)- One 



8 



[ 1 1 1 t 1 t 1 1 1 1 

□ -logD(T) 
> - log D est (T) 




O 

! . S 


1/T 



0.0 0.5 1.0 1.5 2.0 2.5 



FIG. 15: Comparison of the inverse diffusion constant, 
1/D(T), with the prediction l/Z) C5t (T) from Eq. r = 
200. 



ends up with much too small apparent activation ener- 
gies. Again, multi-minima correlations turn out to be 
crucial for the characterization of MB depths. 

In principle, the results of Fig. [l4| may slightly 
change if all MBs rather than those with lifetimes larger 
10 5 MD steps were considered. However, our analysis has 
clearly revealed (see, e.g., Fig. |l|) that the depth of the 
trap only mildly varies when comparing MBs with simi- 
lar e MB . Thus inclusion of MBs with smaller values of r 
would not significantly change the values of the apparent 
activation energies (e MB ). 

Finally, we show that these results, in conjunction 
with p(e MB ;T), largely explain the behavior of the dif- 
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fusion constant D(T). This is a conceptually impor- 
tant step, since we link D(T) to purely structural and 
thermodynamical quantitites, see Eq. 0. The key is the 
mean lifetime (t(c mb ]T)) of MBs at energy e MB , which is 
parametrized by ro(e MB ) and -E ap p(e M B) (Eq. 0). The for- 
mer, to(£ M b), however, has not been deduced from PEL 
properties. Its variation with MB energy is not strong 
(Fig. |2|(c)), so we can hope that setting it to a constant 
will be a good approximation. Thus, Eq. |^ becomes 

D(T) « / d £MB p( £MB ;T) e -^p>MB) = D (T)< 

6Nt J 

(21) 

The estimated diffusion constant derived from this ex- 
pression is shown in Fig. [IJ|. The agreement of D{T) 
with our estimate is satisfactory below T = 1, albeit we 
find a slightly too strong temperature dependence for the 
lowest T. The deviation at T — 2 is due to the depart of 
(T(e MB ; T)) from Arrhenius behavior, see Fig. @(a). 



VI. BARRIER CROSSING. 

When making use of Eq. |l3|, we presumed that the bar- 
riers V(Cab) — y(£,a) in fact are the determinants of the 
temperature dependence of rates. The excellent agree- 
ment between E, pp (e MB ), determined from dynamics, and 
the E° st (t MB ), from the analysis of the PEL, strongly 
indicates that this presumption is indeed true. We will 
show here in a very detailed way that at T = 0.5 = 1.1T C , 
escapes from stable MBs are perfectly activated. More 
precisely, two conditions are fulfilled, (i) the potential 
barriers are much larger than fcsT, (ii) rates follow from 
the ID energy profile of the RP plus corrections from 
perpendicular curvatures. 

We will check these conditions explicitly here by an 
analysis of escape dynamics out of MBs. We made the 
observation that during every escape from a stable MB, 
at least one single barrier larger than 6k^T must be sur- 
mounted. Moreover, this larger jump is mostly under- 
taken from one of the lowest minima of the MB, comp are 
Fig. H From the repeated escape runs of section III B , we 
selected the most frequent ten transitions of that kind. 
From the respective TSs, Q, we computed the RPs, de- 
noted 0( s ); ' = 1-..10. We then investigated the mo- 
tion within the MBs over a long period of the simulation 
where no escape had happened (10 7 MD steps each MB). 
The goal was to observe how the system tries to climb 
the different RPs. To this end, we projected the instanta- 
neous configuration x(t) onto each of the RPs, according 
to 



Sl (t) = {s': \\x(t) - OOOII = mm\\x(t) - Q(s)\\} , 

which means the point on the RP next to x(t). Due to the 
long residences in the MBs, motion therein is largely equi- 
librated. Hence, if the potential energy profiles V ((i{si)) 
along the reaction paths are of importance for the tran- 
sition rates, we expect that the populations pi(s{) of the 




FIG. 16: Parametric plot showing the correspondence of 
— lnp'(s;) + const, to the free energy profile Fi(si)/IcbT + 
const., I = 1...10, T = 0.5. All curves were shifted to start in 
the origin. Insets: comparison of the free energy profiles of 
two reaction paths with the population along the path. 



RPs follow from Boltzmann's law 

pi(si) (xexp{-pV l (s l )}Y l L (s l ) = ex.p{-PF l (s l )}. 

The vibrations perpendicular to path I are accounted for 
by the harmonic partition function 



J dx exp j-^^A^ j<5( 



x ■ i(si)), 



where A„ are the eigenvalues of the Hessian matrix, 
H(si), x v the components of x along the eigenvectors, 
and t(si) is the tangent to the reaction path. 

The upper inset of Fig. |l6| shows an example of pi{si) 
vs. Fi(si) I k-^T . The population of the reaction path 
follows nicely the prediction from its energy profile. For 
RPs with complicated shapes, this correspondence can be 
disturbed. The worst agreement of the considered RPs 
is shown in the second inset. Still, a clear correlation of 
RP population with energy is present. We compiled the 
results for all 10 RPs in Fig. [r] as a parametric plot of 
— In pi vs. Fi/k^T . Curves of slope one result from a per- 
fect equivalence of pi to Fi/k^T . Here, we find an average 
slope of 0.92. Since transition rates are proportional to 
the population of TSs, the implications of these results 
are obvious: MB jump rates follow from energy barriers. 
We finally note that the vibrational terms h\Y l L (si) are 
minor as compared to f3V(s{). 

In view of these results, it is a little surprising that 
the TS location with the help of the auxiliary potential 
V was that unsuccessful (cf. section IV). Since the RP 
population suits well the harmonic description of the RP, 
one expects that motion near the TS is quite harmonic, 
too. Minimizing V in a harmonic potential directly yields 
the stationary state. Consequently, one should easily find 
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the TS when starting from a configuration at s ps 0. 
After section [V, this is not the case, so at least minor 



anharmonicities must be present. 



VII. DISCUSSION. 

The metabasin concept is at the heart of the present 
study. The important insight is that, upon cooling, not 
only the time scale of inter- basin transitions becomes well 
separated from intra- basin vibrations, but also that a 
similar separation occurs between MB hopping and intra- 
MB transitions. Recently, Biroli and Kurchan have anal- 
ysed the general problem of defining metastable states in 
glassy systems jl6[ . They conclude that one has no abso- 
lute notion of a state without making reference to a time 
scale and hence to dynamics. Also the present definition 
of MBs relies on the dynamics of the system. It is, how- 
ever, independent of time scale and exclusively depends 
on the Pback values which directly reflect the topological 
properties of the PEL. 

Our MB definition (see Fig. |l0| ) is devised to eliminate 
the information on trivial back-and-forth jumps within 
MBs. This strongly correlated type of motion is remi- 
niscent of the particles' rattling in the cages formed by 
their neighbors. Similarly, escaping from MBs seems to 
be equivalent to the breaking of cages and thus to struc- 
tural relaxation. Guided by this idea, we have examined 
MB relaxation in great detail: 

First, for repeated relaxation from the same MB, we 
calculated the mean relaxation time (rj) and found Ar- 
rhenius behavior in all cases. The simplest view is 
that the apparent activation energies 2£ app (i) from the 
Arrhenius-like (r,) (Fig. |l|) should correspond to the 
depths of these MBs, i.e. to the typical heights of barriers 
that surround the MBs. Indeed, this has been quantita- 
tively confirmed for the four randomly selected, low-lying 
MBs (see Fig. [ll]). A direct conclusion from the con- 
stancy of E app (i) is that the system does not find smaller 
and ever smaller barriers upon decreasing T. 

Although not of statistical relevance for the whole 
PEL, the results for the four single MBs give us a de- 
tailed picture of the local PEL topography. An important 
outcome is the variation of barrier heights with temper- 
ature, see Fig. [l2| We have already discussed that low 
barriers increase upon cooling, due to enhanced multi- 
minima correlations (growing MBs), while unnecessarily 
high barriers are suppressed. Both effects seem to cancel, 
so that the mean barrier, £"° st , remains constant, lead- 
ing to Arrhenius behavior below T = 1. This cancellation 
may be fortunate, at least we can offer no explanation for 
it, here. As depicted in Fig. [l^, the distribution of bar- 
riers becomes more and more narrow when going from 
T = 0.8 to T = 0.5, but the mean value, i.e. E£ p (l), 
remains constant. The constant apparent activation en- 
ergy of MB 1 down to T = 0.45 implies that the mean 
value of the distribution of barriers has not increased. We 
thus speculate that the growth of barriers due to increas- 



ing multi-minima correlations has essentially come to an 
end at T < 0.5. Although the temperature dependence 
of the barrier distribution has only been analyzed for a 
single MB, the constancy of apparent activation energies 
of the other three MBs and the temperature indepen- 
dence of E app (e M B ) support this idea. Stated differently, 
the development of superstructures of minima seems to 
cease at some temperature above T g . Expressed by p back , 
this means that no minimum with p Wk < 50% will sur- 
pass Pback = 50% upon further cooling, thus being un- 
able to join the MB in question. Hence, an escape se- 
quence found at one temperature T < 0.5 has the same 
length at another one, i.e., from some temperature on, 
the minimum £ m rr) remains at p back < 50% for T — > 0; 
we then say it terminates the sequence. It is an interest- 
ing question under what circumstances such termination 
happens. A trivial example would be a 'transit' minimum 
with one backward and one forward exit, where taking 
the forward leads to a minimum with p back « 0. If the 
backward barrier was higher than the forward one, p back 
would go to zero for T — > 0. On the other hand, the min- 
ima inside MBs generally feature growing p back 's upon 
cooling, because the energetic gain of returning becomes 
more and more attractive. Ideally, thus, for T — > 0, we 
would have p back — » 1 within MBs, and p back — > outside. 
This provides a plausible, physical basis for computing 
barrier heights according to Eq. at least in the limit 
T — » 0. Clearly, a more detailed investigation of the tem- 
perature dependence of p back is necessary to back these 
conclusions. 

Second, we analyzed the average relaxation times 
{t(c mb ;T)) from MBs at fixed energy e MB . Again, they 
displayed Arrhenius behavior, with apparent activation 
energy E app (e MB ) (see Fig. ||), which compared well with 
the prediction from PEL barriers (Fig. |l4). In this con- 
nection, a recent paper by Grigera et aL[B0l is of inter- 
est. The authors use the T^-potential to compute saddles 
in a binary soft-sphere mixture (N — 70). From the 
TSs among these saddles (index one, no shoulder), they 
perform steepest descents to obtain the connected min- 
ima. They define barriers as the energy difference AU 
from the TSs to the lower one of the connected min- 
ima e = min(eo,ei). Plotting the average AC/Ye), they 
find a similar curve to our E app (e MB ), Fig. |l4j, i.e., a 
strong increase of barriers towards lower energies. In 
contrast, when carrying out the same analysis for our 
BMLJ65, we found a nearly constant AU(e), a curve close 
to the first barriers of escapes Eqi shown in Fig. |lj. We 
would have expected this result, since the multi-step na- 
ture of escapes in the BMLJ65 has clearly been demon- 
strated. On the other hand, the contrasting result of 
Grigera et al. indicates that the soft-sphere PEL is not 
organized in multi-minima superstructures. A clarifica- 
tion of this point would be very useful. 

Note that -E app (e M B) is of special importance since it 
bridges the separation between dynamics (diffusion con- 
stant D(T)) and thermodynamics (population of e MB ). 
Clearly, an understanding of -E app (e M B) from basic prin- 
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ciples is highly desirable. It is plausible that the simple 
form E app (e MB ) = e th — e MB can only be expected for a sys- 
tem acting as a completely correlated entity. In contrast, 
two independently relaxing subsystems should generally 
show a weaker dependence of -E app (e MB ) on e MB . This 
can be seen by a very simple argument. Consider 
two independent, identical systems, with MB energies 
£MB, e MB and activation energy E app (e < i^ B ) ). What can 
be said about E app (e MB ) of the union of these systems, 
at MB energy 6 MB — ^Lb ~i~ ^Lb? In the limit of low 
temperatures, the apparent activation energy is given by 
min(i? app (eiJ B ), S app (e^ B )). A proper average over differ- 
ent realizations = e MB - e^l yields S app (e MB ) of 
the combined system. Instead of carrying out this av- 
erage, we use the fact that -E app (e M B) is a monotonous 
function and estimate 

< mm(i4 P ( e «)X»,p(e£ B )) < £C PP (e MB /2). 

Thus, < -E, pp (e MB ) < £ app (e MB /2), which means that 
the combined system shows a weaker dependence on e MB 
than a single copy. For a reasonable PEL topology, one 
would expect |di? app (e M B)/de M B| < 1, because barriers 
should not mount up more than one descents in the 
PEL. Since the e M B-dependence of -E app becomes weaker 
for larger systems, it in turn should increases towards 
smaller N . As a speculation, this might open a way of 
estimating the size of cooperative regions. 

The results shown in Fig. [15], obtained via Eqs. [|and[7], 
demonstrate the use of the present work. From PEL 
barriers (-E app (eMB)) and thermodynamics (p(e MB ;T)) we 
are able to produce a reasonable estimate of dynamics. 
An overall proportionality factor 1 /tq remains as an ad- 
justable parameter, since it could not be predicted from 



PEL structure. As discussed in section [II, one may use 
p(e;T) instead of p(e mB ; T) , since they are nearly identi- 
cal. This is very convenient, because upon contructing 
p(e;T), no information about dynamics is needed. The 
breakdown of the Arrhcnius form of (r(e MB ;T)) above 
2T C limits our description to the temperatures T < 2T C . 
In any event, we would not have dared to make quantita- 
tive statements on the basis of the hopping picture above 
the landscape-influenced temperature regime. At T = 1, 
for instance, we have (e MB ) T ~ —289 (Fig. ||), where 
we already find PEL barriers £ app (-289) « 1 = k B T 
(Fig. |). 

From the fact that we could quantitatively relate MB 
lifetimes to PEL barriers below 2T r and the results from 



section VI , we see that there exist activated barrier cross- 
ing events significantly above T c . As shown before [fl9|| , 
these escape processes from stable MBs determine the 
temperature dependence of the diffusion constant also 
above T c . Thus, the general statement that hopping 
events are not relevant there (see, e.g. 0) is not cor- 
rect for the BMLJ system. This implies that the ideal 
MCT can be applied to systems for which activated pro- 
cesses determine the time-scale of relaxation. Thus it 
seems that the theoretical description of the cage effect 
in terms of structural quantities, as done in MCT, works 



independent of whether the cage effect is purely entropic 
(like in hard-sphere systems) or is to a large degree based 
on barrier-crossing events. 

Moreover, with the help of the unbiased quantity 
p BB (T), we were able to measure the population of basin 
borders. No indication for an abrupt change of relaxation 
mechanism could be observed in p BB (T); in contrast, the 
separation of intra- and inter-basin motion seems to hap- 
pen rather smoothly (see Fig. Q). Thus, there is no qual- 
itative change around T c . 

We finally discuss the relation of our work to the in- 
stantenous normal mode approach (INM) which consid- 
ers the number of 'diffusive modes', / diff , to be at the 
physical basis of diffusion |l3|, |l4|, |l^|. From the direc- 
tions corresponding to negative eigenvalues of the Hes- 
sian H(x{t)) (unstable directions), one filters out the 'dif- 
fusive' directions. Considering the energy profile on the 
straight lines along the unstable directions, La Nave et al. 
observed extremely small barriers, indicating completely 
'entropic' dynamics at the considered temperatures Jl8| . 
This conclusion, though reached for a model of super- 
cooled water, is in contrast to our findings of the rele- 
vance of energetic barriers. A possible key to this ap- 
parent contradiction is that f ditt is directly related to the 
fraction of time spent in 'mobile' regions of configuration 
space. In contrast, we have concentrated on the dura- 
tions of the stable, immobile structures. As the conse- 
quence of longer and longer residences in deep MBs, the 
mobile fraction becomes smaller and smaller. Thus, one 
observes a relation between D(T) and / diff , although it 
is the long trapping times which are the reason for the 
slowing down of dynamics. 

We further note that the MB concept is not imple- 
mented in the INM approach. Supercooled water, e.g., 
exibits very pronounced MB correlations in the time se- 
ries of minima, even for a 'large' system of 216 parti- 
cles p8|| . Generally, fragile glass formers are expected 
to have a 'rugged' PEL, i.e. exibit extensive superstruc- 
tures of minima |2(J. In view of this insight, the success 
of INM analyses for the latter type of systems is quite 
surprising. 



VIII. CONCLUSION. 

Our goal in this paper was to shed some light on the 
temperature dependence of the diffusion constant. In 
our previous work |l{J, metabasins turned out as a use- 
ful concept that reduces correlations of subsequent PEL- 
hopping events. Taking seriously these correlations, the 
present investigation went a step further into this direc- 
tion, by relating the temperature dependence of relax- 
ation to the depths of these multi-minima superstruc- 
tures. We have shown in this paper that a quantita- 
tive link between PEL structure and dynamics is possible 
above T c . However, our approach is still phenomenolog- 
ical at this stage: we are far from predicting -E app (e M B) 
from more general PEL properties or even the interaction 
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potentials. To achieve this is a major challenge, implying 
a deep understanding of energy landscape topology. 

Further insights might be obtained by unveiling the 
real-space aspects of MB relaxation. Here, the corre- 
spondance of MBs to the cage effect should serve as the 
guiding principle. An interesting question along this line 
would be if some of the real-space phenomena found in 



supercooled liquids (e.g. the string-like motion discov- 
ered by Donati et al. |M) can be traced back to energy 
landscape features. 

We are pleased to thank M. Fuchs, D.R. Reichman, and 
H.W. Spiess for helpful discussions. This work has been 
supported by the DFG, Sonderforschungsbereich 262. 
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